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Abstract 

We present a strong lens analysis of SDSS J1004+4112, a unique quasar lens produced by a massive 
cluster of galaxies at z = 0.68, using a newly developed software for gravitational lensing. We find that our 
parametric mass model well reproduces all observations including the positions of quasar images as well 
as those of multiply imaged galaxies with measured spectroscopic redshifts, time delays between quasar 
images, and the positions of faint central images. The predicted large total magnification oi ^^70 suggests 
that the lens system is indeed a useful site for studying the fine structure of a distant quasar and its host 
galaxy. The dark halo component is found to be unimodal centered on the brightest cluster galaxy and 
the Chandra X-ray surface brightness profile. In addition, the orientation of the halo component is quite 
consistent with those of the brightest cluster galaxy and member galaxy distribution, implying that the 
lensing cluster is a relaxed system. The radial profile of the best-fit mass model is in good agreement with a 
mass profile inferred from the X-ray observation. While the inner radial slope of the dark halo component 
is consistent with being —1, a clear dependence of the predicted A-D time delay on the slope indicates 
that an additional time delay measurement will improve constraints on the mass model. 

Key words: dark matter — galaxies: clusters: general — galaxies: quasars: individual (SDSS 
J1004-I-4112) — gravitational lensing 
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1. Introduction 

SDSS J1004-I-4112 is a unique quasar lens system (Inada 
et al. 2003; Oguri ct al. 2004; Inada et al. 2005; Inada ct al. 
2008). A quasar at z = 1.734 is multiply imaged into five 
images, with their maximum image separation of 14'.' 7, 
produced by a massive cluster of galaxies ai z = 0.68. It 
is one of only two known examples of cluster-scale quasar 
lenses, the other being the triple lens SDSS J1029-K2623 
with the maximum image separation of 22."5 (Inada et al. 
2006; Oguri et al. 2008). In addition to the quasar images, 
SDSS J1004-I-4112 contains spectroscopically confirmed 
multiply imaged galaxies at z ^ 3 (Sharon et al. 2005). 
Both the quasar images and the lensing cluster have been 
detected in Chandra X-ray observations (Ota et al. 2006). 
Moreover, time delays between some of the quasar images 
have also been measured from long-term optical monitor- 
ing observations (Fohlmeistcr ct al. 2007; Fohlmeister ct 
al. 2008). 

Such a wealth of observational data available enable de- 
tailed investigations of the central mass distribution of the 
lensing cluster. Indeed, there have been several attempts 
to model the mass distribution of SDSS J1004-I-4112, us- 
ing either parametric or non-parametric method. Oguri 
et al. (2004) adopted two-component (halo plus central 
galaxy) model to successfully reproduce the positions of 
four quasar images, but even the parities and temporal 
ordering of the quasar images could not be determined 
because of model degeneracies. Kawano & Oguri (2006) 



extended mass modeling along this line, and explored how 
time delay measurements can distinguish different mass 
profiles. Fohlmeister et al. (2007) pointed out that it is 
important to include perturbations from member galaxies 
to reproduce the observed time delay between quasar im- 
ages A and B. On the other hand, Williams & Saha (2004) 
and Saha et al. (2007) performed non-parametric mass 
modeling to show possible substructures in the lensing 
cluster. From the non-parametric mass modeling, Saha 
et al. (2006) and Liesenborgs et al. (2009) concluded that 
the radial mass profile is consistent with that predicted in 
A''-body simulation (Navarro et al. 1997). We summarize 
previous mass modeling in Table 1. 

In this paper, we revisit strong lens modeling of 
SDSS J1004-f-4112 adopting a parametric mass model. We 
include many observational constrains currently available, 
including central images of lensed quasars and galaxies, 
fiux rations, and time delay between quasar images (see 
Table 1). In particular this paper represents first para- 
metric mass modeling that includes both the quasar time 
delays and the positions of multiply imaged galaxies as 
constraints. We compare our best-fit mass model with 
the Chandra X-ray observation of this system (Ota et al. 
2006). 

This paper is organized as follow. We describe our mass 
model in §2. We show our results in §3, and give conclu- 
sion in §2. A new software for gravitational lensing, which 
is used for the mass modeling, is presented in Appendix. 
Throughout the paper we adopt the matter density of 
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Table 1. Previous mass modeling of SDSS J1004+4112 



[Vol. , 



Reference 



model* 



constraints (quasar) ^ constraints (galaxy)' 



Inada et al. (2003) 
Oguri et al. (2004) 
Williams & Saha (2004) 
Sharon et al. (2005) 
Kawano & Oguri (2006) 
Saha et al. (2006) 
Fohlmeister et al. (2007) 
Saha et al. (2007) 
Inada et al. (2008) 
Liesenborgs et al. (2009) 



SIE+pcrt 


pos+fiux (4) 




NFW+SIE+pert 


pos-f flux (4) 




non-parametric 


pos (4) 




SPL+gals 


pos (5) 


pos (5) 


gNFW+SIE+pert 


pos+fiux (4) 




non-parametric 


pos (5) 


pos (8) 


NFW+deV+gals+pcrt 


pos-(-flux (5) 




non-parametric 


pos (5), At (2) 


pos (8) 


NFW+SIE+gals-t-pert 


pos+fiux (5), At (2) 




non-parametric 


pos (5), At (3) 


pos (7) 



gNFW+Jaffe+gals+pcrt pos+fiux (5), At (3) 



pos (27) 



This work 

* Name of models: SIE = singular isothermal ellipsoid, pert = external perturbation (e.g., external shear), NFW = 
Navarro, Frenk & White (NFW) profile, SPL = softened power-law model, gals = perturbations from member galaxies, 
typieally modeled by truncated isothermal profiles, gNFW = generalized NFW profile, dev = de Vaucouleurs profile, Jaffe 
pseudo-Jaffe profile. 



t 

parentheses show the number of images used as constraints. 



'pos" indicates constraints from image positions, "flux" from fluxes of images, and At from time delays. Numbers in 



^M — 0.26 and the cosmological constant of rt\ — 0.74. 
but regard the dimensionless Hubble constant ft, as a pa- 
rameter. With this choice of cosmological parameters, a 
physical transverse distance oilh^^ kpc at the redshift of 
the lensing cluster (2 = 0.68) corresponds to 0.20 arcsec. 
We denote a angular diameter distance from observer to 
lens as Di, from observer to source as Dg, and from lens 
to source as Dig. 

2. Mass Modeling 

2.1. A Parametric Model 

We model a main halo of the lensing cluster as the gen- 
eralized NFW profile (e.g., Jing, Suto 2000). Its radial 
density profile is given by 

Ps 



p{r) 



(1) 



(r/r,)"(I + r/r,)3-"' 

In this model, the inner slope is parametrized by a (0 < 
q; < 2) ; the original NFW profile corresponds to a = I . We 
adopt tire following modified concentration parameter as 
a model parameter: 



C-2 = 



r-2 



(2 - a)r. 



(2) 



where r_2 indicates the radius where the radial slope be- 
comes dlnp/dlnr = —2 and rvu- is the virial radius of the 
cluster. The characteristic density is described as 

/\{z)p{z)c' 



3mg„fw(c) 

TOgnfw(c) 



(1 



-dr. 



(3) 



(4) 



where A(z) is nonlinear overdensity at redshift z which 
we adopt values predicted by the spherical collapse model. 
The lensing deflection angle is related with the projected 
two-dimensional mass distribution (i.e., convergence k) 
computed by 



n{r) 



1 



?( V r^ 



z'^)dz, 



(5) 



^crit 

with Ecrit = {c'^/4:TtG){Ds/D\D\s) being the critical surface 
mass density computed from angular diameter distances 
between observer, lens, and source. 

We then introduce an ellipticity in the projected mass 
distribution by replacing r in K(r) by the following quan- 
tity: 



K(r) : 



(1-e) 



+ (l-e)y2 



(6) 



where e is an ellipticity (the axis ratio is 1 — e), and x and 
y are defined by 

X = xcos9e+ysin6e, (7) 

y = — xsin6'e + 2/cos6'e. (8) 

In this paper we take the a;-axis to West and the y-axis 
to North, and therefore the position angle 9e is measured 
East of North. 

We include the brightest cluster galaxy (BCG) Gl mod- 
eled by pseudo-Jaffe Ellipsoid (Keeton 2001b). The con- 
vergence takes the following form: 



K = Kq 



1 



1 



y^s'^/q + x'^ +y'^/q'^ Y'^^V^+'^^'+^V^ 



kq : 



,(9) 



(10) 



2GE„,tA'V9' 

where q ~ 1 — e is the axis ratio, s is the core radius, a 
is the truncation radius, and a is the velocity dispersion. 
This profile has a constant density at r ^ s, an isothermal 
density profile r^^ at s<C?'^a, and a sharply decline 
profile r~^ at r ^ a. 

Perturbations from member galaxies are also included. 
We model individual galaxies by the pseudo-Jaffe ellipsoid 
without core radius (eq. [9] with s = 0). We fix positions, 
relative luminosities, ellipticities and position angles of 
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Table 2. Constraints from Icnsod quasar images 



Table 3. Constraints from Icnsed galaxy images 



Name 


Ax n 


Ay n 


Am 


At [days] 


A 


0.000 


0.000 


= 


= 


B 


-1.317 


3.532 


0.35±0.3 


-40.6 ±1.8 


C 


11.039 


-4.492 


0.87±0.3 


-821.6±2.1 


D 


8.399 


9.707 


1.50±0.3 




E 


7.197 


4.603 


6.30 ±0.8 





I'he quasar redshitt is Zs = 1.734. i'hc positional error is assumed 
to 0."04 for all the Icnscd quasar images. 

member galaxies to observed values, and assumes that 
the velocity dispersions and truncation radii scale with 
the luminosity: 

1/4 

(11) 



a 
a. 



1/2 



(12) 



The mass-to-light ratio becomes constant with this scal- 
ing. We include 14 member galaxies within ^ 20" from 
the center, which are selected from the gri-hand Subaru 
Suprime-cam images (Oguri et al. 2004). We adopt r-band 
luminosities for the scaling. 

To achieve better fit, we also include additional several 
perturbations. We consider general perturbations whose 
lens potentials cfi are described as (see, e.g., Evans & Witt 
2003; Kawano et al. 2004; Congdon & Keeton 2005; Yoo 
et al. 2006) 



r'^cosm{9 -9^- 7r/2). 



(13) 



In this paper we include four perturbation terms with m = 
2 (external shear; e.g., Keeton et al. 1997), 3, 4, and 5. 

2.2. Observational Constraints 

We adopt positions of five quasar images measured by 
Inada et al. (2005) using the Hubble Space Telescope 
Advanced Camera for Surveys (HST/ACS) F814W im- 
age. Considering possible effects of microlensing or small- 
scale structure, we adopt conservative positional errors of 
0''04, and also relative magnitudes of 0.3 (0.8) for image 
B-D (E). In addition, we include measured time delays 
between image A and B (Fohlmeister et al. 2007) and be- 
tween image A and C (Fohlmeister et al. 2008). When 
fitting the time delays, we allow the Hubble constant to 
vary with a Gaussian prior of /i = 0.72 ± 0.04. 

We also include multiply imaged galaxies, identified by 
Sharon et al. (2005), as constraints. We revisit deep multi- 
band HST/ACS images (F435W, F555W, and F814W; 
GO-10509 and GO-10716), and identify several features 
associated to each lensed images. We use positions of all 
these features for our mass modeling. We include central 
images as well (Liesenborgs et al. 2009). We assume larger 
positional errors of 0''4 than those of the quasar images, 
partly because the determination of the centroids of the 
extended galaxy images are much less accurate. 

Figure 1 shows the HST/ACS image of 



Name 



Aa- 



Ay 



Al.l 
A1.2 
A1.3 
A1.4 
A1.5 
A2.1 
A2.2 
A2.3 
A2.4 
A2.5 
A3.1 
A3.2 
A3.3 
A3. 4 
A3.5 
Bl.l 
B1.2 
B1.3 
B2.1 
B2.2 
B2.3 
Cl.l 
C1.2 
C1.3 
C2.1 
C2.2 
C2.3 



3.33 



3.33 



3.33 



2.73 



2.73 



3.28 



3.28 



3.93 


-2.78 


1.33 


19.37 


19.23 


14.67 


18.83 


15.87 


6.83 


3.22 


4.13 


-2.68 


1.93 


19.87 


19.43 


14.02 


18.33 


15.72 


6.83 


3.12 


4.33 


-1.98 


2.73 


20.37 


19.95 


13.04 


18.03 


15.87 


6.83 


3.02 


8.88 


-2.16 


-5.45 


15.84 


8.33 


2.57 


8.45 


-2.26 


-5.07 


16.04 


8.33 


2.57 


10.25 


-3.06 


-7.55 


15.39 


8.49 


2.72 


9.95 


-3.36 


-7.30 


15.44 


8.49 


2.72 



The positional error is assumed to 0."4 
for all the lensed galaxy images. The 
redshifts are measured spectroscopi- 
cally. 



SDSS J1004+4112, together with the positions of 
multiple images summarized in Tables 2 and 3. A notable 
feature of this cluster strong lens system, which can 
easily be seen in the Figure, is that multiple images are 
distributed in a very wide range in radius, ranging from 
the central images very near the cluster center to lensed 
galaxy images at ^ 30" from the cluster center. This is 
apparently good for constraining the density profile of 
the lensing cluster. 

We also add several Gaussian priors to the mass model. 
Based on the measurement by Inada et al. (2008), we as- 
sume the velocity dispersion of the central galaxy Gl to 
a = 352 ± 13 kms^^.^ The position of Gl is fixed to the 
observed position, (7"114, 4"409) in our coordinate sys- 
tem whose origin is at the quasar image A. From the ob- 
served position and shape, we assume the ellipticity and 
the position angle to 0.3 ±0.05 and 152 ±5 deg, respec- 
tively. Furthermore we add a weak prior to the truncation 



Strictly speaking, the velocity dispersion cTobs computed from 
the density profile can in principle differ from the input param- 
eter (T for the pseudo Jaffc profile. However, from Eliasdottir et 
al. (2007) we find that a ^ ctq^s for values similar to those in the 
best-fit model {s/a Ri 0.05), which suggests that our assumption 
of (T = fT|-,i3g is reasonable. 
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Fig. 1. Left: The HST/ACS image of SDSS J1004+4112. North is up and West is right. Stellar objects labeled by A-E are the 
5 lensed quasar images. A large galaxy superposed on image E is the brightest cluster galaxy Gl. Right: Positions of the 5 lensed 
quasar images {filled circles) and lensed galaxies images (other symbols). The size and orientation of the panel is same as the left 
panel. Different symbols have different rcdshifts. See Tables 2 and 3 for the relative coordinate values. 



radius, a = 



' ± 4", based on the observed correlation be- 



tween the velocity dispersion and the truncation radius 
(Natarajan et al. 2009). 

2.3. Model Optimization 

We use a software named glafic for all the calculations of 
lens properties and model optimizations (see Appendix 1). 
We employ a standard x^ minimization to find the best- 
fit mass model. Specifically we adopt a downhill simplex 
method to find a minimum. To speed up the calculations, 
we estimate x^ in the source plane, which is found to be 
sufficiently accurate for our purpose. See Appendix 2 for 
detailed discussions about the source plane x^ minimiza- 
tion. 



3. Result 

3.1. B est- fit NFW Model 

First, we fix the inner slope of the dark halo component 
to a = 1 (i.e., the NFW profile) and derive the best-fit 
mass model. Figure 2 shows best-fit critical curves. It 
is seen that the best-fit model reproduces the observed 
multiple image families quite well. In addition, it success- 
fully reproduces observed time delays between quasar im- 
ages. The best-fit model has x^ = 31 for 39 degree of free- 
dom, suggesting that our choice of errors was reasonable. 
The contribution from each source is reasonably similar, 
X^ = 4.7 from the quasar, 21 from the galaxy A, 0.9 from 
galaxy B, and 2.6 from galaxy C. The best-fit model pre- 
dicts magnifications of the five quasar images to ha ~ 29.7, 
Hb = 19.6; ^iQ = 11.6, /xd = 5.8, and /xe = 0.16. The to- 
tal magnification for all the quasar images is /Xtot = 67. 
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Fig. 2. The best-fit mass model assuming the NFW {a = 1) 
profile for the dark halo component. Symbols are same as 
those in the right panel of Figure 1, but for showing the image 
positions predicted by the best-fit model. Thick solid lines in- 
dicate critical curves for the quasar source redshift, Zs = 1.734, 
whereas thin lines are critical curves for the highest redshift 
among the multiply imaged galaxies, Zs = 3.33. 
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The model also predicts the time delay between quasar 
image A and D to be A^ad = A^d — At a = 1218 days, 
and that between quasar image A and E to be A^ae = 
AiE - AiA = 1674 days. The AD time delay is slightly 
smaller than the lower limit reported by Fohlmeister et 
al. (2008), At AT) > 1250 days. 

We find that the best-fit centroid of the dark halo 
(NFW) component is (6.92^0.32, 4.25to:'^4) at 95% con- 
fidence limit, which is quite consistent with the observed 
position of Gl. The result is in marked contrast with 
Oguri et al. (2004), in which significant offsets between 
the center of the dark halo component and that of Gl 
have been reported based on modeling of quasar images 
A-D. Such large offset is no longer allowed because of ad- 
ditional constraints from multiply imaged galaxies, time 
delays, and central images. The best-fit center of the dark 
halo component is also consistent with the X-ray center 
reported by Ota et al. (2006). Furthermore, the best-fit 
position angle of 9^ = 152.9 deg for the dark halo compo- 
nent is quite consistent with the position angle of Gl, and 
also that of the member galaxy distributions studied in 
Oguri et al. (2004). The concentricity and alignment be- 
tween dark matter, BCG, and X-ray implies that the lens- 
ing cluster is a highly relaxed system (see also Liesenborgs 
et al. 2009). The best- fit ellipticity of the halo component 
is e = 0.26. 

The best-fit parameters for perturbations terms are (e, 
6/,)=(0.040, 51.8) for m = 2, (0.019, 114) for m = 3, (0.013, 
47) for m = 4, and (0.010, 16.5) for m = 5. Thus the per- 
turbations are rather small, but they are still important 
for accurate reproduction of lensed images, particularly 
for those of the quasar (see also Oguri et al. 2004). 

One of the most important quantity to characterize 
strong lensing system is the Einstein radius rEin- We com- 
pute the Einstein radii rsin for our best-fit mass model 
using the following relation: 



M{< rEi, 



''Ein^crit- 



(14) 



We find TEin = 8."14 for the quasar redshift Zg = 1.734, and 
13"38 for the redshift of the lensed galaxy A, Zg = 3.33. 
If we compute rEin only from the dark halo compo- 
nent excluding any contributions from galaxies, we obtain 
rEin = 4."84 for Zs = 1.734 and 10"31 for Zg = 3.33, which 
are quite different from those compute from the total mass 
distribution. This suggests that the effect of the BCG Gl 
on the lens system is quite significant. 

3.2. Comparison with X-ray 

Next we compare the best-fit radial mass profile derived 
from strong lensing with that inferred from the Chandra 
X-ray observation (Ota et al. 2006). In brief, from the 
Chandra observation the extended X-ray emission from 
the lensing cluster was detected out to ~ lf5 from the clus- 
ter center, with the temperature of ^ 6.4 keV. Assuming 
the isothermal profiles, Ota et al. (2006) constrained the 
projected mass profile and argued that the mass within 
100 kpc agrees well with the mass expected from strong 
lensing. Here we compare our result of new mass modeling 
with the X-ray result. 
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Fig. 3. Projected two-dimensional cumulative mass distri- 
butions from lensing and X-ray. Thick and thin lines are 
best-fit total and dark matter mass distributions from strong 
lens mass modeling, respectively. Open circles arc mass dis- 
tributions inferred from X-ray surface brightness and temper- 
ature measurements by Chandra, assuming an isothermal gas, 
/3 model for the gas profile, and hydrostatic equilibrium (Ota 
et al. 2006). Errors are from the temperature measurement 
uncertainty. 

Figure 3 compares the projected two-dimensional mass 
profiles from gravitational lensing and X-ray measure- 
ments. We confirm that the profiles agree quite well with 
each other, including radial slopes of the profiles. The 
agreement suggests that the effect of the halo triaxiality, 
which affects the apparent two-dimensional lensing masses 
particularly near the center of the cluster (see, e.g., Oguri 
et al. 2005; Gavazzi 2005), is not significant. 

However, it should be noted that the best-fit halo 
mass of Mvir = 1.0 x 10^^h~^MQ and the concentration 
of c_2 = 2.8 are quite different from those inferred from 
X-ray, Mvh- - 4.3 x lO^^/i-^f© and c„2 ^ 6.1. One rea- 
son for this the strong degeneracy between Afvir and c_2 
inherent to strong lens mass modeling. Basically strong 
lenses constrain the central core mass of the cluster, which 
is a strong function of both A/vir and c_2. The determi- 
nation of Mvir and c_2 solely from strong lensing relies on 
the extrapolation of the subtle change of the radial slope 
out to much larger radii. The robust determination of 
these parameters from lensing therefore requires addition 
constraints from weak gravitational lensing. 

3.3. Generalized NFW Profile 

Now we allow the inner slope of the dark halo compo- 
nent a to vary to see how well the current strong lens data, 
including time delays which are quite helpful to break de- 
generacy in mass models (e.g., Kawano & Oguri 2006), 
can constrain the inner density profile. Specifically, for 
each fixed value of a we optimize the other parameters. 
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Fig. 4. Upper: The x difference, Ax =x X (iniii): ^^ ^ 
function of the inner slope of the dark halo component. For 
each a, the other model parameters are optimized to mini- 
mize x^. The solid line shows the default case, whereas the 
dashed line indicates the result when an additional prior of 
C—2 = 6.0 ±1.5 from the Chandra X-ray observation is in- 
cluded for mass modeling. Lower: The total magnification 
factor for the five quasar images, /itoti and the time delay 
between the quasar images A and D, AJaDi predicted by the 
best-fit models as a function of the inner dark halo slope a. 
The horizontal dotted line indicates the observed lower limit 
of AiAD reported by Fohlmeister et al. (2008). 



Figure 4 shows the x^ difference as a function of a. We 
find that our mass modchng is quite consistent with the 
NFW profile, i.e., a = 1. Wc constrain the range of the 
slope to 0.76 < a < 1.41 at 95% confidence limit. The 
profile as steep as a = 1.5 is clearly rejected. There is a 
clear correlation between a and the velocity dispersion of 
galaxy Gl such that the best-fit velocity dispersion de- 
creases with increasing a, which approximately conserves 
the central core mass of the total matter distribution. 

As discussed in § 3.2, the Chandra X-ray observation 
suggests that the lensing cluster may have larger value of 
the concentration parameter than the best-fit NFW model 
predicts. We include this effect by adding a prior of c_2 = 
6 ± 1.5 to the mass model to see how the constraint on a is 
modified. The result shown in Figure 4 indicates that the 
constraint is basically shifted to the lower a, i.e., shallower 
inner density slope. The resulting range is 0.62 <a < 1.14 
at 95% confidence limit. 

In Figure 4, wc also show the total magnification fac- 
tor for the five quasar images, /Xtot, and the time delay 
between the quasar images A and D, AiAD, predicted by 
the best-fit model for each fixed a. We find that /itot 
is decreasing and AiAD is increasing with increasing a, 
which are consistent with well-known dependence of the 
magnification and time delay on the radial density slope 
(e.g., Wambsganss & Paczynski 1994; Oguri & Kawano 
2003). Thus the reported lower limit of A^ad > 1250 days 
(Fohlmeister et al. 2008) prefers a steeper inner slope 
(larger a), which appears to be opposed to the effect of 
the prior from X-ray discussed above. In either case, the 
strong dependence of A^ad on a implies that additional 
A-D time delay measurements will greatly help to con- 
strain the mass model further. 

4. Summary 

We have revisited parametric mass modeling of 
SDSS J1004-I-4112, a unique quasar-cluster lens system, 
using a newly developed mass modeling software. We in- 
clude several new constraints, including positions of spec- 
troscopically confirmed multiply imaged galaxies, time de- 
lays between some quasar images, and faint central im- 
ages. Our model comprising of a halo component mod- 
eled by the generalized NFW profile, member galaxies 
including the brightest cluster galaxy Gl, and perturba- 
tion terms, well successfully reproduced all observations 
including time delays. 

Unlike earlier claims based on parametric mass model- 
ing, we have found that the center and orientation of the 
dark halo component is in good agreement with those of 
member galaxies and Chandra X-ray observation, imply- 
ing that the cluster is highly relaxed. The radial profile 
from strong lensing is also in excellent agreement with 
the mass profile inferred from the X-ray observation. Our 
mass modeling prefers the dark halo component with the 
inner slope close to a — 1, being consistent with so-called 
NFW density profile. Additional measurement of the time 
delay between quasar image A and D will be useful to 
constrain the mass model further. The predicted total 
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magnification of /xtot = 67 for the NFW profile is quite 
large compared with those for typical galaxy-scale lenses, 
because of the shallower density profiles for clusters. This 
makes the lens system a good site for studying the fine 
structure of the quasar through microlensing (Richards et 
al. 2004; Green 2006; Lamer et al. 2006; Pooley et al. 2007) 
or for studying host galaxies of distant quasars (Ross et al. 
2009). We note that our result based on parametric mass 
modeling is broadly consistent with recent non-parametric 
mass modeling by e.g., Liesenborgs et al. (2009). 

Our result suggests that the core of the lensing clus- 
ter at z = 0.68 is highly evolved. Recently, Limousin et 
al. (2010) showed that the cluster MACSJ1423.8-h2404 at 
z = 0.54 is similarly highly relaxed based on the compar- 
ison of mass, light, and gas distributions. Therefore our 
result may point to the fact that relaxed clusters are quite 
common already at z '--^ 0.6. 

It is clear that the lensing cluster SDSS J1004+4112 
is currently one of the best-studied high-rcdshift clusters 
whose inner density profile is very tightly constrained by 
strong lensing. Additional constraints on this cluster with 
weak lensing, Sunyaev-Zel'dovich effect, and spectroscopic 
identifications of member galaxies will be important to 
advance our understanding of high-rcdshift clusters. 

I would like to thank the referee, Marceau Limousin, 
for useful comments and suggestion. 

Appendix 1. Lens Modeling Software glafic 

We have developed a comprehensive software pack- 
age called glafic that can be used for a wide va- 
riety of analysis for gravitational lensing. Its fea- 
tures include efficient computations of Icnscd images 
for both point and extended sources, handling of mul- 
tiple sources, a wide range of lens potentials, and 
the implementation of noble technique for mass mod- 
eling. Currently the software can be downloaded from 
http: //www. slac . Stamford. edu/~oguri/glaf ic/. 

In this code, we adopt the adaptive mesh refinement 
algorithm described in Kceton (2001a) for deriving lensed 
point source images, although the code uses rectangular 
coordinates rather than polar coordinates. Critical curves 
are computed using a marching squares technique (Julio 
et al. 2007). In simulations of extended sources, one can 
convolve point spread functions and include a number of 
galaxies read from a catalog file, which should also be 
useful for weak lensing work. For more details, readers 
are referred to a manual available at the URL above. 

Appendix 2. Source-plane x^ Minimization 

Strong lens modeling with the standard x^ minimiza- 
tion is sometimes time-consuming, especially when many 
lens potential components and images are involved. One 
way to overcome this problem is to evaluate x^ ^^ the 
source plane instead of the image plane. Although the 
source plane x^ involves approximations (given that ob- 
servational measurements are always made in the image 



plane) and therefore is less accurate than the image plane 
X^, it allows one to estimate x^ without solving the non- 
linear lens equation. This technique has been adopted 
by several authors (e.g., Kayscr 1990; Kochanck 1991; 
Koopmans et al. 1998; Keeton 2001a; Smith et al. 2005; 
Bradac et al. 2005; Halkola et al. 2006; Julio et al. 2007; 
Sand et al. 2008), although the implementations were 
quite different for different papers. Here we describe our 
implementation and argue the accuracy of the technique. 
For a given source position u, x^ is estimated as 
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where x^, ^i, Ati is the position, magnification, and time 
delay for the i-th image, respectively. They are related 
with the source position through the lens equation: 



V0(x, 
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Here 4>{x.) is the so-called lens potential. From the lens 
equation, the magnification factor /i; is computed as 



^i = |dct(M,)|, 
with Mi being the magnification tensor defined by 

1 du 
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Finally the time delay Ati is 
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Computing x^ in equation (Al) usually requires the 
derivation of x^ for a given u, which is the most time- 
consuming part because the lens equation is not one-to- 
one mapping and thus the extensive solution finding in the 
image plane is needed. In the source plane x^ technique, 
we do not solve the lens equation but just approximate 
X^ as follows. Assuming Xj and x,;.obs are close with each 
other, we can write 

X,,obs - Xi W M, (Ui^obs - u) , (A9) 

where u^^obs is the source position computed from the ob- 
served i-th image position: 

U,,obs = ^i,ohs - V(?!)(Xi_obs), (AlO) 

and Mi is estimated at x = Xi^obs- Then equation (A2) 
becomes 

(u,,obs - u)'^Mf (Ui^obs - U) 
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Similarly, the magnification and time delay are approxi- 
mated as 
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(relative) time delay error of 0.5 days. We run a Markov 
Chain Monte Carlo (MCMC) for this model to explore 
X^ around the best-fit model, and at each step we com- 
pute both Xhng ^iid Xsrc to See the difference of these. 
The result shown in Figure 5 indicates that Xsrc a-grees 
well with XhiiK 'within a few percent level, which arc suf- 
ficiently accurate to derive the best-fit mass model and 
^'^|iv3;*.\^ errors on best-fit parameters. We note that Xg^^, is even 
more accurate when observational constraints are tighter. 
We note that Halkola ct al. (2006) also studied accuracy 
of the source plane x^ based on modeling of the strong lens 
cluster A1689. They found a clear correlation between 
X?jj and Xsrc t)^t with different values. This is because 
they did not take account of the full magnification tensor 
(eq. [A7]) but simply adopted a magnification factor to 
approximate x^i as has been often done in the literature. 
1 O Thus our result highlights the quantitative importance of 

the proper approximation taking the full lens mapping 
into account. 
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Fig. 5. Accuracy of the source plane x , Xsrc- Each point 
shows the fractional difference between xfmE '^'^'^ Xsrc ^.s a 
function of v? , estimated at each MCMC step. The as- 
sumed mass model is described in the text. 
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In the code, d^/dx is evaluated numerically. Inserting 
these expressions into equations (A3) and (A4) yield 
Xflux.src ^^'^ XAt.src respectively. The source plane x^ 
is just a sum of these three: 

(A14) 
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Note that ttiq and Atg are nuisance parameters whose 
best-fit values can easily be derived as 
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Also note that one can adopt image fluxes as constraints 
rather than magnitudes. The modification for this is quite 
straightforward (see also Keeton 2001a). 

We find that the source plane x^j if properly evalu- 
ated like above, is accurate enough to be used in most 
cases of strong lens mass modeling. As a specific exam- 
ple, we consider a mass model consisting of an SIE and 
external shear. The model parameters are a = 320kms^^, 
e = 0.3, ^e = 20 deg, 7 = 0.1, 6*^ = 90 deg, zi = 0.5, z, = 2.0, 
u=(0.04, —0.02). We then add errors to the image po- 
sitions, magnitudes, time delays, and fit the "observed" 
four-image system using the same mass model. We adopt 
a positional error of 0."01, a magnitude error of 0.1, and a 
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